#ifndef _EFEM_SF_H
#define _EFEM_SF_H

#include "exception"
#include "stdexcept"
#include "cmath"
// Include shape functions of scalar values
#include "nfem_sf.h"

#ifndef EFEM_SF_ZERO
#define EFEM_SF_ZERO 1e-10
#endif // EFEM_SF_ZERO

// Space of edge based finite element method vector shape (basic) function
namespace efem_sf
{

enum valtype_e
{
    Value, // Original vector components
    Gradient, // Gradients of the vector components
    Curl, // Curl components
};

enum ordertype_e
{
    Normal, // indice of the edge's nodes are ordered from small to big
    Reverse, // indice of the edge's nodes are ordered from big to small
};

/**
 * @brief      The 1st order (linear) vector basis functions
 */
class linear
{
public:
    linear(){}
    virtual ~linear(){}

    /**
     * @brief       Edge based vector shape function defined on a 2D triangle.
     * 
     * This function has three nodes on the corners, respectively.
	 * Indice of the nodes and edges are as show.
	 * 
	 * y
	 * |
	 * |
	 * 2
	 * | \
	 * |   \
	 * |    \
	 * |      \
	 * |       \
	 * |         \
	 * 0-------------1--->x
     *   edge0: 0->1
     *   edge1: 1->2
     *   edge2: 2->0
     * 
     * @param[in]  x      Evaluation position x
     * @param[in]  y      Evaluation position y
     * @param[in]  x1     x coordinate of the first node of the triangle
	 * @param[in]  x2     x coordinate of the second node of the triangle
	 * @param[in]  x3     x coordinate of the third node of the triangle
	 * @param[in]  y1     y coordinate of the first node of the triangle
	 * @param[in]  y2     y coordinate of the second node of the triangle
	 * @param[in]  y3     y coordinate of the third node of the triangle
	 * @param[in]  idx    Edge index
     * @param[in]  vt     Calculation target. Could be Value, Gradient or Curl
     * @param[in]  ot     Ordering direction of the edge's nodes
     * @param      xval   Evaluated shape function. Equals to Vx (Value), dVx/dy (Gradient) or curl(V) (Curl)
     * @param      yval   Evaluated shape function. Equals to Vy (Value), dVy/dx (Gradient) or curl(V) (Curl)
     */
    void triangle(double x, double y, double x1, double x2, double x3, 
		double y1, double y2, double y3, unsigned int idx, valtype_e vt, 
        ordertype_e ot, double &xval, double &yval);

    /**
	 * @brief    Edge based vector shape function defined on a 3D tetrahedron.
	 * 
	 * This function has four nodes on the corners, respectively.
	 * Indice of the nodes and edges are as show.
     * 
	 * q
	 * |
	 * 3
	 * |\
	 * |  \
	 * |	\
	 * |	  \
	 * |		\
	 * |		  \
	 * 0------------2---> v
	 * |         /
	 * |       /
	 * |      /
	 * |    /
	 * |   /
	 * | /
	 * |/
	 * 1
	 * |
	 * u
     * 
     * edge0: 0->1
     * edge1: 0->2
     * edge2: 0->3
     * edge3: 1->2
     * edge4: 1->3 (or 3->1 in some other implementations, Does not really matter ^,^)
     * edge5: 2->3
	 *
	 * @param[in]  x      Evaluation position x
	 * @param[in]  y      Evaluation position y
	 * @param[in]  z      Evaluation position z
	 * @param[in]  x1     x coordinate of the first node of the tetrahedron
	 * @param[in]  x2     x coordinate of the second node of the tetrahedron
	 * @param[in]  x3     x coordinate of the third node of the tetrahedron
	 * @param[in]  x4     x coordinate of the fourth node of the tetrahedron
	 * @param[in]  y1     y coordinate of the first node of the tetrahedron
	 * @param[in]  y2     y coordinate of the second node of the tetrahedron
	 * @param[in]  y3     y coordinate of the third node of the tetrahedron
	 * @param[in]  y4     y coordinate of the fourth node of the tetrahedron
	 * @param[in]  z1     z coordinate of the first node of the tetrahedron
	 * @param[in]  z2     z coordinate of the second node of the tetrahedron
	 * @param[in]  z3     z coordinate of the third node of the tetrahedron
	 * @param[in]  z4     z coordinate of the fourth node of the tetrahedron
	 * @param[in]  idx    Edge index
	 * @param[in]  vt     Calculation target. Could be Value, Curl or Div
     * @param[in]  ot     Ordering direction of the edge's nodes
     * @param      xval   Evaluated shape function. Equals to Vx (Value), dVx/dy (Gradient) or curl(V)_x (Curl)
     * @param      yval   Evaluated shape function. Equals to Vy (Value), dVy/dx (Gradient) or curl(V)_y (Curl)
     * @param      zval   Evaluated shape function. Equals to Vz (Value), dVz/dx (Gradient) or curl(V)_z (Curl)
     * @param      xval2  Evaluated shape function. Equals to dVx/dz (Gradient)
     * @param      yval2  Evaluated shape function. Equals to dVy/dz (Gradient)
     * @param      zval2  Evaluated shape function. Equals to dVz/dy (Gradient)
	 */
    void tetrahedron(double x, double y, double z, double x1, double x2, double x3, double x4, 
		double y1, double y2, double y3, double y4, double z1, double z2, double z3, double z4,  
		unsigned int idx, valtype_e vt, ordertype_e ot, double &xval, double &yval, double &zval, 
        double &xval2, double &yval2, double &zval2);
};

/**
 * @brief      The 2nd order (quadratic) vector basis functions
 */
class quadratic
{
public:
	quadratic(){}
	virtual ~quadratic(){}

	/**
     * @brief       Edge based vector shape function defined on a 2D triangle.
     * 
     * This function has three nodes on the corners, respectively.
	 * Indice of the nodes and edges are as show. The shape functions 
	 * are defined on the edges and the interior facet
	 * 
	 * y
	 * |
	 * |
	 * 2
	 * |\
	 * |  \
	 * |    \
	 * |      \
	 * N4       N3
	 * |          \
	 * |            \
	 * |      N7       \
	 * N5     |-> N6     N2
	 * |                   \
	 * |                     \
	 * 0-----N0->-----N1->-----1--->x
     *   edge0: 0->1
     *   edge1: 1->2
     *   edge2: 2->0
     * 
     * @param[in]  x      Evaluation position x
     * @param[in]  y      Evaluation position y
     * @param[in]  x1     x coordinate of the first node of the triangle
	 * @param[in]  x2     x coordinate of the second node of the triangle
	 * @param[in]  x3     x coordinate of the third node of the triangle
	 * @param[in]  y1     y coordinate of the first node of the triangle
	 * @param[in]  y2     y coordinate of the second node of the triangle
	 * @param[in]  y3     y coordinate of the third node of the triangle
	 * @param[in]  idx    function index (0->7)
     * @param[in]  vt     Calculation target. Could be Value, Gradient or Curl
     * @param[in]  ot     Ordering direction of the edge's nodes. Could be Normal or Reverse
     * @param      xval   Evaluated shape function. Equals to Vx (Value), dVx/dy (Gradient) or curl(V) (Curl)
     * @param      yval   Evaluated shape function. Equals to Vy (Value), dVy/dx (Gradient) or curl(V) (Curl)
     */
    void triangle(double x, double y, double x1, double x2, double x3, 
		double y1, double y2, double y3, unsigned int idx, valtype_e vt, 
        ordertype_e ot, double &xval, double &yval);
};

}

#endif // _EFEM_SF_H